function [price,outRate,inRate]=snowball_MC(S0,C,knockinB,knockoutB,r,Sigma,T,Ntraj,Nstep)
%%%%%%%%%%%%%函数说明%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%传统蒙特卡罗方法计算雪球产品价格%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%输入参数：
%%%%%%%%%%S0：标的初始价格；knockinB：敲入界限；knockoutB:敲出界限；C:年化票息率；
%%%%%%%%%%r：无风险利率；Sigma：标的波动率；T：到期时间；
%%%%%%%%%%Ntraj：模拟路径数；Nstep：每条路径一年的时间分割数
%%%%%%%%%%输出参数：
%%%%%%%%%%Price:雪球产品价格;outRate:敲出路径所占比例;inRate:敲入未敲出路径所占比例
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EXAMPLE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<1
T=1;S0=100; r=0.03; Sigma=0.13; knockinB=0.8;knockoutB=1.03;C=0.2;
Ntraj=100; Nstep=200;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=simuPath(S0,T,r,Sigma,Ntraj,Nstep);%一列为一条价格路径
payoff=zeros(1,Ntraj);
holdingPeriod=T*ones(1,Ntraj);%记录雪球产品的有效持有期

%Case I: knock out
outORnot=S>S0*knockoutB;%检验节点价格是否出现敲出
outPath=find(any(outORnot,1));
period=[0,1:T*Nstep]'.*outORnot(:,outPath);%用于寻找敲出路径的有效持有期
period(find(period==0))=NaN;
holdingPeriod(outPath)=min(period,[],'omitnan')/Nstep;
notoutPath=setdiff([1:Ntraj],outPath);%未敲出
% Case II: knock in but not out
inORnot=S<=S0*knockinB; % filter knock in
inPath=find(any(inORnot,1));
innotoutPath=setdiff(inPath,intersect(inPath,outPath));% filter in but not out
% case III: filter not in not out
notinnotoutPath=setdiff([1:Ntraj],[inPath,outPath]);%未敲入未敲出

%%%%%%%计算每一种情况的收益
payoff(outPath)=S0*(1+C*holdingPeriod(outPath));
payoff(notinnotoutPath)=S0*(1+C*holdingPeriod(notinnotoutPath));
payoff(innotoutPath)= S0*min(S(Nstep,innotoutPath)./S0,1);

price=exp(-r*T)*mean(payoff);%计算所有路径收益均值并进行贴现
outRate=length(outPath)/Ntraj;
inRate=length(inPath)/Ntraj;

% visualization
disp(['nums of out path is ', num2str(length(outPath))]);
disp(['nums of knock in but not out path is ', num2str(length(innotoutPath))])
disp(['nums of not in not out path is ', num2str(length(notinnotoutPath))])
end